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We study nonequilibrium steady states of lattice gases with nearest-neighbor interactions that 
are driven between two reservoirs. Density profiles in these systems exhibit oscillations close to the 
reservoirs. We demonstrate that an approach based on time-dependent density functional theory 
copes with these oscillations and predicts phase diagrams of bulk densities to a good approximation 
under arbitrary boundary-reservoir couplings. The minimum or maximum current principles can be 
applied only for specific bulk-adapted couplings. We show that they generally fail to give the correct 
topology of phase diagrams but can still be useful for getting insight into the mutual arrangement 
of different phases. 
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Lattice gases provide useful models for investigat- 
ing driven particle transport in biological, chemical and 
physical systems [J, H L They are used to study Brownian 
'3, organic photovoltaic cells @, 
6] , to name only some recent ap- 
plications. In these examples, the asymmetric simple ex- 
clusion process (ASEP) appears as a basic building block 
for the description of low-dimensional transport. As such 
it has developed into one of the standard models for in- 
vestigating nonequilibrium steady states (NESS). In the 
ASEP, particles hop under the influence of a bias field and 
cannot occupy the same place. For open boundary condi- 
tions, these biased systems show intriguing phenomena, 
like boundary- induced phase transitions 7], emergence 
of shock fronts and self-organized pattern formation 
. To explain these phenomena, it is generally sufficient 
to consider the totally asymmetric simple exclusion pro- 
cesses (TASEPs), where the particle transport in the bias 
direction is unidirectional. 

Much progress has been made in the past to under- 
stand ASEPs and TASEPs (for reviews, see Refs. [13- 
12[), but only a few studies so far have addressed TASEPs 
with particle-particle interaction going beyond hard-core 
repulsions Blil-El- While in this case a complete de- 



scription of the NESS seems to be out of reach for open 
systems, it has been shown that the minimum/maximum 
current principles [?[ or the domain- wall (shock- 
front) theory [8|, [ij, [iSj , can be used in a specific model 
setup to predict phase diagrams of the bulk density. 
This setup requires that the relations between correla- 
tion functions and densities close to the system-reservoir 
boundaries are the same as in the bulk. To ensure this, 
a particular way of particle injection and ejection has to 
be taken, which was used in Refs. l44l7|. 

However, realistic reservoirs are not of that type but 
are specified by only a few control parameters, as, for ex- 
ample, temperature and chemical potential. It is there- 
fore needed to develop methods that can deal with such 
situations both on basic reasons and in view of various 



applications. For example, in driven transport through 
molecular bridges (e.g.. in the incoherent limit for weak 
coupling to leads) [19|, |20| , one does not specify a com- 
plicated bath-system coupling, but is led by the fact that 
the typical relaxation dynamics in the bath is much faster 
than in the system. Accordingly, the baths are supposed 
to be in equilibrium and characterized essentially by their 
chemical potential. 

In this Letter we consider a TASEP with nearest- 
neighbor interactions as sketched in Fig. [1] and show 
that the couplings of the system to the reservoirs have a 
strong influence on the density profiles in the nonequilib- 
rium steady state (NESS). In general, oscillations of these 
profiles occur close to the boundaries, similar as they are 
known for equilibrium systems. We present a theory that 
is able to cope with these oscillations and to describe the 
density profiles to a good approximation. Based on the 
theory, the phase diagram of the bulk density can be de- 
termined. Our theoretical approach provides a general 
method to predict such phase diagrams for driven sys- 
tems with interactions under general boundary-reservoir 
couplings. 

To demonstrate our theoretical approach we consider 
as an example the following specific model, see Fig. [T] 
Particles with repulsive nearest-neighbor interaction of 
strength V > perform a unidirectional hopping motion 
between neighboring sites of a lattice with N sites and 
are thereby transported from a left to a right particle 
reservoir. The microstate of the system is specified by 
the set of occupation numbers n = {n^}, i = 1, . . . , N, 
where = or 1 if the corresponding site is vacant 
or occupied by a particle. The jump rate of a particle 
from a site z to a vacant neighboring site (i + 1) is F^ = 
exp(— A£^/2), where AE is the energy difference (in units 
of the thermal energy) between the final and initial state 
after and before the jump. 

When V > Vc = 2.89, the bulk current-density relation 
of this model exhibits a double-hump structure with two 
maxima at densities Pmax,i and Pma.x,2, and a minimum 



in between at p = 0.5 17|. The bulk dynamics are con- 
sidered to apply to jumps from sites i = 2, . . . , (N — 2), 
while couplings to the reservoirs are taken into account 
for injections from and ejections to the reservoirs. In ad- 
dition we need to specify rates for jumps from sites i = 1 
and {N — 1), where, if one adopted the bulk rates, a 
nearest- and a next-nearest-neighbor site would be miss- 
ing to the left and right, respectively. As sketched in 
Fig. [11 two cases of boundary couplings are considered. 
These couplings are referred to as "bulk- adapted" and 
"equilibrated-bath" coupling and explained for the left 
reservoir in the following. Corresponding couplings are 
applied to the right reservoir. 

For the equilibrated-bath coupling [Fig. [rfb)], the left 
reservoir is considered to be an equilibrated ideal Fermi 
gas with a chemical potential /Xl, corresponding to a 
reservoir density Pl = l/[exp(— /Xl) -|- 1]. Accordingly, 
we write for the rate Fl of particle injection V 1^(112) = 
PLexp[(/XL-ri,2y)/2], and Fi(n3) = exp(-n3F/2) for the 
jump rate from site one. 

The bulk-adapted coupling [Fig. [Ila)] is arranged in 
such a way that the system can be viewed as being 
continued into a reservoir with density Pl, correspond- 
ing to relations between correlation functions (riinj . . .), 
and densities pi — (n^) as in the bulk, where (...) de- 
notes an average over the distribution of microstates 
in the NESS. In a closed bulk (ring) system, when an 
initial configuration {ni+i = 0,ni+2} would be given, 
two rates are possible for a particle jump from site i 
(i.e. Ui = 1): Fi = exp(-ni+2y/2), if 0, 
while Fi — exp[(l — ni+2)W2]i if "-i-i = 1- For given 
{ni_|_i = 0,ni_|_2}, let us denote by p(01|0?T.i+2) and 
p(ll|0ni+2) the conditional probabilities for the config- 
urations {rii_i,ni} = {0,1} and {rii_i,ni} — {1,1} 
to occur in the NESS of a closed bulk system with 
density p^ and interaction V ^ respectively. The in- 
jection rate FL(n2) then results from a weighting of 
rates with the probabilities p(01|0n2) and p(ll|0n2) cor- 
responding to virtual configurations {n_i = 0, no = 
l,ni = 0,71-2} and {n_i = l,r7,o = l,ni = 0,ri2} at 
the boundaries, i.e., FL(n2) = p(01[0n2) exp(— 712^/2) -|- 
^((ll|07^2)exp[(l - n2)Vl2\. Analogously, Fi(n3) = 
p(0[10n3)exp(-?i3F/2)-Hp(l[10n3)exp[(l-n3)V^/2]. For 
a practical implementation of a corresponding kinetic 
Monte Carlo (KMC) simulation, the conditional prob- 
abilities p(.|.) are determined from separate KMC sim- 
ulations of periodic ring systems with a particle density 
/9 = Pl as indicated in Fig. [iTa) . The rates for the two 
models are summarized in [2l|. 

KMC simulations have been carried out for both the 
bulk-adapted and the equilibrated-bath coupling and re- 
sults are shown by the squares and circles in Fig. [2] 
for V — 2Vc, Pl — 0.9 and pn — 0.7. As a conse- 
quence of the specific arrangements in the bulk-adapted 
case, monotonously varying density profiles are obtained, 
similar as in the TASEP with hard-core repulsion only. 
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FIG. 1. Sketch of the model and illustration of the particle 
injection for (a) the bulk-adapted, and (b) the equilibrated- 
bath couplings. Right reservoirs are not shown. 



The value of the bulk density pe = 0.70 agrees with 
that predicted from applying the maximum current prin- 
ciple to the bulk current-density relation. For the 
equilibrated-bath coupling by contrast, pronounced oscil- 
lations appear at the boundaries. This implies that meth- 
ods relying on the bulk current-density relation (mini- 
mum/maximum current principles, shock- front/domain- 
wall theories) cannot be applied any more to predict 
the bulk densities. In fact. Fig. [5] shows that a value 
Pb = 0.40 is obtained for the equilibrated-bath coupling, 
which differs from that for the bulk-adapted coupling. 

The question is whether a theory can successfully ac- 
count for the bulk densities and associated phase dia- 
grams for interacting driven particle systems that gen- 
erally will show density oscillations at the boundaries. 
Mean-field theories with simple factorization schemes, as, 
for example, {niHij) ~ {ni){n.j) fail for the model intro- 
duced here, since they are not even capable to predict 
the double-hump structure in the bulk current-density 
relation. We now show that the application of the time- 
dependent density functional theory (TDFT) presented 
in [lli nil can deal with the complications associated 
with the density oscillations. In this approach the av- 
erage current ji from site i to site (i -I- 1) in the bulk 
(i = 2, . . . , (A^ - 2)) is given by 



{Pi+2 - Ci+i)e + pi+i - pi+2 + C'i+i 



Pi - Cj 
PiPi+i 



(1) 



where pi — {nt) and p^ = 1 — (rii) = 1 — pi are the particle 
and hole density, respectively, and Ci = {nirii+i) is the 
two-point correlation function [23j , 

a = exp i-V) — -— ;- . (2 

1 - Pi - Pj+1 + Ci 

This expression can be explicitly solved to yield functions 
Ci = Ci{pi^ pi+i) and in this way, the bulk currents ji in 
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FIG. 2. (Color online) Steady-state density profiles for V — 
2Vc from KMC simulations (symbols) in comparison with the 
TDFT predictions (lines). For the KMC results averages were 
performed over 10^ particle jumps in the steady state. The 
inset shows the density profile for the equilibrated-bath cou- 
plings close to the left and right boundary, respectively. 



Eq. ([T]) become functionals of the density profile p{t) = 

For tlie bulk-adapted couplings, Eq. ([T]) applies also at 
the boundaries, that means for i = 0, 1, {N — 1) and 
N 24|. For the equilibrated-bath couplings, however, 
different functional structures are obtained: 



Jo = Phe 



JN-I 



Pi-C, 
h 

Pn-1 — Cn-1 
Pn-1 



ip3 - C2) + P2-P3 + C'2 



JN = PrS 



Cn-i£^^^ + Pn — Cw-1 



Pn~i — 



(3a) 
(3b) 

! 1 

(3c) 
(3d) 



To compare the predictions of this theory with the 
KMC results we have integrated numerically the cou- 
pled set of rate equations dpi{t)/dt = ji-i{t) — ji{t) 
{i = 1, ...,iV) with the currents given by ([I|)-(I21) and 
evaluated the long-time limit to analyze the NESS. As 
shown in Fig. [21 the density profile from the TDFT for 
the bulk-adapted coupling agrees well with KMC results. 
Surprisingly, also the oscillations of the KMC density 
profiles for the equilibrated-bath couplings are closely re- 
produced by the TDFT, see the inset of Fig. [2] Moreover, 
the TDFT gives a value pa = 0.42 that is only slightly 
larger than the corresponding KMC value = 0.40. 
For the steady-state currents we find Jness = 0.174 in 



the KMC and j^ess = 0.149 in the TDFT. 

With the density profiles determined from either KMC 
or TDFT we can identify the singularities in the depen- 
dence of the bulk density on the reservoir densities Pl 
and Pr. Corresponding transition lines are shown in the 
phase diagram of Fig.|3]for V — 2T4, and for (a) the bulk- 
adapted and (b) the equilibrated-bath couplings. Over- 
all, the TDFT accounts well for the phase transitions as 
determined from the KMC simulations. Due to the con- 
struction of the bulk-adapted couplings, the knowledge 
of the exact bulk current-density relation in the NESS 
would allow one to determine exactly the phase diagram 
by applying the minimum/maximum current principles. 
In this way seven phases are identified in Fig.[3l^a). Since 
the TDFT does not yield the bulk current-density rela- 
tion in the NESS exactly, small deviations are seen in 
Fig. [Sla) between KMC and TDFT resuhs. 

The strongly different diagram in Fig. ^h) shows that 
the minimum/maximum current principles are no longer 
successful. By contrast, the TDFT accounts well for the 
five phases identified in the KMC simulations. There 
are three phases, where ps is some function of either Pl 
or Pr, and a maximum current phase with pe = Pmax,i 
and a minimum current phase with pe — 0.5. Generally, 
Pb must be either determined by Pl or p^, or by the ex- 
trema in the bulk current-density relation, since the mini- 
mum/maximum current principles still apply in the inner 
bulk regions, where the density profiles are monotonously 
varying. The absence of a maximum current phase with 
Pb = Pmax,2 for the equilibrated-bath couplings means 
that also the topology of the phase diagram in Fig. [Sjb) 
is changed compared to that in Fig.[3l[a). It can be shown 
that this change of topology is generally possible and that 
there are no further hidden phases in addition to the ones 
shown in Fig. E^b) |25j]. Those phases, which appear in 
Fig-Etb), have a connection to the ones in Fig.[3l[a) in the 
sense that their mutual arrangement remains the same. 

In conclusion we have shown that the theoretical ap- 
proach based on the TDFT can cope with the problem of 
driven lattice gases with extended interactions, where os- 
cillations in densities and correlation functions naturally 
occur at the boundaries. The oscillations imply that the 
minimum/maximum current principles are no longer suf- 
ficient, since they can only be applied in an inner bulk 
region, where the density profile varies monotonously. 
When, for unmodified bulk dynamics, deliberately chang- 
ing the boundary couplings to the bulk-adapted ones 
in order to enable the use of the minimum/maximum 
current principles, the currents and bulk densities in 
the NESS as well as the associated phase diagrams are 
strongly infiuenced. The principles are nevertheless use- 
ful to justify the identification of possible phases on the 
basis of a known bulk current-density relation. The in- 
vestigation of the bulk-adapted couplings can be helpful 
to predict the mutual arrangement of those phases that 
appear for the boundary couplings of interest. 
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FIG. 3. Phase diagrams of the NESS for (a) the bulk-adapted 
and (b) the equihbrated-bath coupHngs a.t V = 2Vc. KMC 
results for the phase transitions are marked by the symbols 
and TDFT results by the lines (solid lines for first-order and 
dashed lines for second-order phase transitions). In (a) seven 
phases are obtained, where the bulk density pe equals either 
the reservoir densities Pl or Pr, or the value 0.5 in the mini- 
mum current phase, or the two possible values pmax,i or pmax,2 
in the maximum current phases. In (b) a maximum cur- 
rent phase with pe ~ Pmax,i, a minimum current phase with 
Pb ~ 0.5, and three phases with reservoir-controlled densities 
Pb = /i(Pl), Pb = /iii(Pr), and ps = /v(pr) are obtained, 
where the different functions of Pl or Pr are determined from 
the simulated or calculated density profiles. The labeling by 
roman numbers in (a) and (b) has been chosen in such a man- 
ner that corresponding phases have equal numbers. 



We believe that our theory can be useful to explore a 
wider range of driven open systems which are of impor- 
tance for many biological processes and electronic trans- 
port phenomena in small molecular devices that can be 
treated within the incoherent classical limit. 

We thank W. Dieterich and A. Nitzan for very illumi- 
nating discussions concerning this work. 
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